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Abstract The aim of this paper is to show how option prices in the Jump- diffusion models, mainly on the Merton 
and Kou models, can be computed using meshless methods based on Radial Basis Function (RBF) interpolation. The 
RBF technique is demonstrated by solving the partial integro-differential equation (PIDE) in one-dimension for the 
American vanilla put and the European vanilla call/put options on dividend-paying stocks. The radial basis function we 
select is the Cubic Spline. We also propose a simple numerical algorithm for finding a finite computational range of an 
improper integral term in the PIDE so that the accuracy of approximation of the integral can be improved. Moreover, 
we use a numerical technique called factorization of the Cubic Spline to avoid inverting the ill-conditioned Cubic Spline 
interpolant. Finally, we will show numerically that in the European case the solution is second order accurate for the 
spatial and time variables, while in the American case it is second order accurate for spatial variables and first order 
accurate for time variables. 
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Introduction 



In this paper we show how to compute European and American option prices in the Jump- 
diffusion model using Radial Basis Function (RBF) interpolation techniques. RBF methods 
have recently been proposed for numerically solving initial value and free boundary problems 
for the classical Black and Scholes equation, both in the one and in the multiple asset case 



, 2004a|b 


Hon and Mao 


1999, 


Larsson et al. 


2008) 



present paper is that in the Jump-diffusion model, as in general Levy type models, the Black 
and Scholes PDE is replaced by a Partial Integro-Differential Operator or PIDE, involving a 
global term in the form of an integral operator. The PIDE has a form: 



d T u(x, t) 
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(cf. Cont and Tankov 2004 Schoutens , 2003 ; 2006 ) . Our main contribution is to show how to 



numerically solve ([lj) in an efficient way using RBFs, both for initial value and free boundary 
problems (as for American options). We have chosen the Jump-diffusion model as a typical 
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case on which to test the present RBF methodology. Our method extends however without 
problems to other contexts in which the basic pricing equation is a PIDE, like that of Levy- 
type models such as Carr-Geman-Madan-Yor (CGMY) (Carr et al. 2002 ) or Variance Gamma 
(VG) ( |Carr et al.\ p98j |Madan and Milne[ |l99~I] ) These will be treated in a future paper. 
Currently, PIDEs such as the Merton Model ( |Merton[ |1976[ ) and the Kou Model flKou 



2002[ Kou and Wang, 2001), one have mostly been treated by a traditional Finite Difference 
Method (FDM) or Finite Elements Method (FEM). In FDM, the idea is to simply fully 
discretize the PIDE on an equidistant grid, after having (artificially) localized the equations to 
some bounded interval/domain in M. The global integral term can be computed by numerical 



quadrature or by using the Fast Fourier Transform (FFT) (see, Almendra, 2004; Almendral 



and Oosterlee 


2004} 


2005; 2006; 


2007[ Andersen and Andreasen[ 


2000; 


Briani et al.\ 


2007 


Cont and Volte 


ikova 


2005 


|d'Hal] 


uin et al. 20041 12005[ |Hirsa and Madam 20041 Wang et al. 



2007). By contrast, FEM is defined as piecewise polynomial functions or wavelet functions 



on regular triangularizations. This technique is used to approximate solutions of the partial 
differential terms as well as of the integral term (cf. |Almendral and Oosterlee 2005, Matache 



et al, 2003 2005) 



In general, there is a problem which arises with these current approaches. Some of the 
literature, e.g. (Andersen and Andreasen, 2000 Briani et al., 2007 Cont and Voltchkova 



2005), plays down the importance of pricing American and European vanilla option values 



when time to maturity is less than 3 months. The reason is that for short times- to- maturity 
the numerical methods used to price the option tend to be inaccurate near the strike price 
where a singularity (kink) exists. A singularity is defined as a point at which the function, or 
its derivative, is discontinuous. The payoff functions of vanilla call and put options have such 
a singularity. As a result, standard numerical methods such as FDM with Crank-Nicolson and 
without any adaptive schemes cannot ensure accuracy of option prices around the strike and a 
substantial amount of oscillation occurs around the strike when Option Delta A and Gamma V 



are approximated (Giles and Carter, 2006). Giles and Carter shed light on this kind of problem 



(Giles and Carter 



2006) by suggesting Rannacher's time stepping method. This is a mixture 



of four half-timsteps of backward Euler and Crank-Nicolson methods. Although they solve an 
one dimensional PDE under the Black-Scholes model and Heston's volatility model rather than 
a PIDE under Levy models, their methods of using backwards Euler timestepping in one or 
more initial timesteps have been proved to be achieved second-order convergence in a European 
case. They also carry out a detail error analysis of their methods by using Fourier analysis 
and find out four half-timesteps of backward Euler time-marching is the minimum require to 



recover second-order convergence of solving the PDE. Forysth et ai.(d'Halluin et al. 



also use the similar idea by suggesting Rannacher's time stepping method (Rannacher 



2005) 



1984) 



to solve a PIDE under the Merton Jump-diffusion model. They demonstrate this technique 
by approximating an option price whose maturity is a quarter of a year. This method gives 
second order rates of convergence when pricing European options but not American ones. 
By using the same idea and combining it with a penalty method and a modified form of a 
timestep selector suggested in (Johnson 1987), Forysth et al. (d'Halluin et al. 2004| ) show how 
to achieve second order convergence for pricing American options. Unfortunate they do not 
carry out any stability analysis when they apply Rannacher's time stepping method to solve 
the PIDE. Moreover there is no minimum requirement of choosing half-timsteps of backward 
Euler before Crank-Nicolson methods are applied. All they do is by trial and error. 
In most recent research papers in quantitative finance (cf. Fausshauer et al. 2004a|b Hon 



and Mao, 1999; Larsson et al., 2008) RBF-approximation methods with Multiquadric (MQ) 



as a basis function have been proposed for numerically solving the classical Black and Scholes 
PDE, both in the one and in the multiple asset case. In this literature MQ is a more favorite 
choice than other radial basis functions, such as Thin plate Spline and the like, because of its 
comparatively higher accuracy. MQ contains a shape parameter which plays an imperative 
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role in the accuracy of the method (cf. Wendland , 2005 ) . Most of this recent literature still 
chooses this parameter by trial and error or some other ad-hoc means. Although there exists 
a substantial literature on choosing an "optimal" shape parameter in MQ, e.g. (Fasshauer 



and Zhang, 2007), (Fornberg and Wright, 2004) and (Kansa and Carlson, 1992), it is still an 



open question and there is no theoretical proof for selecting an optimal shape parameter (cf. 
Wendlandf 2005) in MQ. Besides this, the standard approach to the solution of the radial 



basis function interpolation problem has been recognized as an ill-conditioned problem for 
many years (cf. Fasshauer 2007, chapter 16). This is especially true when infinitely smooth 
basic functions such as MQ or Guassian are used with small values of their associated shape 
parameters. More recently, Fasshauer and Mccourt's least-squares approximation based on 



early truncation of the kernel expansion (Fasshauer and Mccourt, 2010 )and Fornberg and 



co-workers's Contour-Pade integration method (e.g. Driscoll and Fornberg, 2002 Fornberg 



and Wright, 2004; Larsson and Fornberg, 2005) are successful in solving the ill-conditioning 



problem of RBF, but the techniques are only restricted to solve the simple interpolation 
problem rather than to solve PDEs, especially parabolic PDEs. Although Ling and his co- 



workers (e.g. Ling and Kansa, 2004; Brown et al, 2005 Ling and Kansa, 2005) address the 



ill-conditioning problem by using preconditioning methods and extend them to solve PDEs, 
the methods are not possible to be applied to solving PIDEs. 

Our RBF-approximation method with the Cubic Spline as a basis function will circum- 
vent these disadvantages. This paper is divided into five sections, including this introduction. 
Section [2] is a brief review of both the Merton and Kou Jump-diffusion models. In section [3] 
we first explain and then define our RBF algorithm for solving PIDEs, which we implement 
the Jump-diffusion model. Section [4] contains our numerical results for both European and 
American call and put options, including an analysis of the max error, the root-mean-square 
error, the rate of convergence and the approximation of A and T and also a comparison the 
accuracy of our solution with that of FDM and FEM . Section [5] concludes. 



2. PIDE Option Pricing Formula in Jump-diffusion Market 

In this short section we will focus on the Merton and the Kou Jump-diffusion Models which 
are general Levy processes consisting of Brownian motion and compound Possion jumps. By 
using these models we can describe the price dynamics of the underlying risky asset, (St)t>o- 
The evolution of (<S*)t>o is driven by a diffusion process, punctuated by jumps which describe 
rare events such as crashes and/or drawdowns at random intervals. As a market model, it is 
an example of an incomplete market. 
The stock price process, (St)t>o, driven by these models, is given by: 

St = S e L < (2) 

where So is the stock price at time zero and L± is defined by: 

N t 

Lf.= j c t + aW t + J2 Y i> ( 3 ) 

i=i 

here, 7 C is a drift term, a is a volatility, Wt is a Brownian motion, Nt is a Possion process 
with intensity A, Yi is an i.i.d. sequence of random variables. Since a > in pi), there exists a 
risk-neutral probability measure Q such that the discounted process {e~^ r ~^>St}t>o becomes 



a martingale (cf. Sato 1999, Theorems 33.1 and 33.2), where r is the interest rate and q is the 



dividend rate. For a discussion of the issue of choosing Q see, for example, ( Cont and Tankov 



2004). Then under this new measure Q, the risk-neutral Levy triplet of Lt can be described 
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as follows: 



(7c cr, v) 



where 



7c 



Xrj + xv(dx), 



(4) 



Here we focus on the case where the Levy measure is associated to the pure-jump component 
and hence the Levy measure v{dx) can be written as Xf(x)dx, where the weight function f(x) 
can take two forms: 

(1) In the classical Merton model, for any i £ {1,2, . . .}, Y{ are log-normally distributed 
variables with Yi ~ N(fXj, a 2 ) and as a result, 



(2) In the Kou model, 



pa\e 



-^ x l x > + (I -p)a 2 e a ^ l x < . 



(5) 



(6) 



Remark 1 : In the Merton Jump-diffusion model, one should notice that Yi is i.i.d so for 
each i S {1, 2, 3, . . .}, Yi has the same mean and variance. For the sake of simplicity, we use 
fij and cjj to represent the mean and variance of each Yi respectively. 

Also in Q, 7] = j R (e x — l)f(x) dx represents the expected relative price change due to a 
jump. Since we have defined the Levy density function f(x) for both Jump-diffusion processes, 
rj can be computed as: 



(1) In the Merton model, 



(2) In the Kou model, 



V 



Hj+(t 2 j/2 



pai (l-p)a 2 



+ 



1. 



a\ — 1 a 2 + 1 

This is found by integrating e x over the real line by setting a\ > 1 and a 2 > 0. 



(7) 



(8) 



For the detai l s of the computation of ([7[) and (|8 b, we shall refer the reader to (Cont and 



Tankov 



2004 



Boyarchenko and Levendorskii 2002). 
The drift-term j c in (|3l) assumes that e~ (r ~ q>t St is a martingale with respect to the natural 



filtration. We let r = T — t, the time-to-maturity, where T is the maturity of the financial 
option under consideration and we introduce x = log St, the underlying asset's log-price. If 
u(x, t) denotes the values of some (American and European) contingent claim on St when 
log = x and r = T — t, then it is well-known, see for example, (Cont and Tankov, 2004) 
that u satisfies the following PIDE in the non-exercise region: 



d T u(x, t) = ^<J 2 d 2 x u + (r — q— ^a 2 - rj } <) r ii - f r + \ )u + 



X / u(x + y,r)f(y)dy, 
Jm 

--: £[u](x,t). 



(9) 
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with initial value 

u(x,0) = g(x) := G(e x ) 



maxje 1 — K, 0} , call option 
maxji^ — e x , 0} , put option 



(10) 



For an American put, we have to take into account the possibility of early exercise (e.g., 



2004; 


Schoutens 


2003; 


2006) 



option can be achieved by maximizing over all allowed exercise strategies: 



U{X,T) 



esssup T , er(tjT) E t Q [e- r ^*- t )G(e x -) 



fill 



where T(t, T) denotes the set of non-anticipating exercise times t* , satisfying t < t* < T. 
To actually compute the u(x, r) of the American put, one can solve the following linear 
complementarity problem (Cont and Tankov, 2004 Schoutens, 20031 2006): 



d T u(r, x) - Cu(x, t) > 0, in (0, T) x E 
u(x, t) - G(e x ) > 0, a.e. in (0, T) x 
(u(x, t) - G{e x )) (d T u(r, x) - £u(x, r)) = 0, in (0, T) x R 

u (x,0) = G(e x ), 



(12) 
(13) 
(14) 
(15) 



Since we only deal with a jump-diffusion model with a > and finite jump intensity in this 
paper, we know that by Pham (Pham, [1997 ), the smooth pasting condition, 

du{x T * , r*) 



dx 



-1 



is valid at time of exercise r*. Therefore the value of an American put option is continuously 
differentiable with respect to the underlying on (0, T) x M; in particular the derivative is 
continuous across the exercise boundary. 

Remark 2 : One should notice that if we set A = 0, Q will become original Black-Scholes 
PDE. 



3. Meshfree Numerical Approximation Method 



Meshfree radial basis function (RBF) interpolation is a well-known technique for reconstruct- 
ing an unknown function from scattered data. It has numerous applications in different fields, 
such as terrain modeling in geology, surface reconstruction in imaging, and the numerical 
solution of partial differential equations in applied mathematics. In particular, RBFs have 
recently been used to solve the PDEs of quantitative finance. A number of authors, including 
Fausshauer et al. ( Fausshauer et al. , 2004a|b ) , Larsson et al. ( Larsson et al. , 2008 ) , Petters- 



son et al. (Pettersson et al., 2008) and Hon and Mao (Hon and Mao, 1999), have suggested 



RBFs as a tool for solving Black-Scholes equations for European as well as American options. 
This numerical scheme for the estimation of partial derivatives using RBFs was originally 



proposed by Kansa (Kansa, 1990a), resulting in a new method for solving partial differential 



equations (Kansa, 1990b). The aim here is to obtain a RBF approximation of the initial value 



or pay-off of the option. Once we are disposition of such an RBF-interpolant, we implement 
an RBF-scheme to solve the PIDE with this RBF-interpolant as initial value. The general 
idea of the proposed numerical scheme is to approximate the unknown function u(x, r) by 
an RBF-interpolant using the interpolation points found for the initial value using the RBF- 
scheme, and derive a system of linear constant coefficient ODE by requiring that the PIDE 
([9]) be satisfied in the chosen RBF-interpolation points. After picking interpolation points 
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Xj E R, we approximate, for any fixed time-to-maturity r, the solution u(x, r) in Q by its 
RBF-interpolant : 

JV 

u(x,t) ~ ^2pj(T)(j)(\\x - Xj\\ 2 ) =: u(x,t), (16) 

3=1 

Since the radial basis function does not depend on time, the time derivative of u(x, r) in 
equation ^ is simply: 

T 3=1 T 

Moreover, the first and second partial derivatives of u(x, r) with respect to x are 

8u(x,r) » d<K\x- Xj \) 
ox z — ' ox 

3=1 

d») 

3=1 

where for the particular case when cp is the Cubic Spline, 

d<p(\x — Xj\) J 3(|x — Xj\) 2 if x — Xj > 0, 



dx 1 — 3(|aj — Xj\) 2 if x — xj < 0, 



d 2 4>(\x — Xj 

dx~ 2 = "VP - ^ 



(20) 



Jiy - 6(|x - xA). (21) 



In this research we choose the Cubic Spline rather than the most popular ones, MQ and IMQ 
as a basis function because of its simplicity and accuracy and without containing any shape 
parameters. 



3.1 Transforming PIDE to A System of ODEs by RBF 



Given a set of interpolation points x±, . . . 



N x N matrices A, A x and A xx defined by 



3 ' ' 



we can construct 
Xj\))^, ,^» r and 



Xi 



, xjv in R, and an RBF 
\Xi - x j\))l<ij<Ni {4> (\xt ^3un<i,j<N 
c j I)) i<i -<n res P ec fi ve ly- Note in case the chosen according to the Equally 

Spacing Method, ESM, used in (Fausshauer et al., 2004a|b Hon and Mao 1999). In brief, 
Equally Spacing Method is the way to choose equally spaced points in a finite interval. In the 
ESM, we determine an interval [x m i n ,x max ] outside of which we can neglect the contribution 
of u(x, t) to the global integral term of a PIDE S, and for given N = 0, 1, 2, ... , simply put 



x 



3 ■- 



x 



Ax 



Xv\ 



+ jAx, j = 0,1,2,..., iV_i 



(22) 



where Ax = (x max — x m \ a )/(N — 1). We also define a matrix- valued function y — > A(y) by 
(<p(\xi + y — x j\)) j <N - If we substitute u(x,t) for u(x,t) in (9) and require the PIDE to 
be satisfied in the interpolation points Xj, we arrive at the following system of ODEs for the 



October 26, 2011 2:28 Applied Mathematical Finance JumpDf'RBFSv4 



7 



vector p(r) : = 



[pi(t),...,pn(t)) 

2/2 

Ap T = y A xxP + ( r - q - y - A?? ) v4 x p + (r + A)Ap + 



A(y)f(y)dy ) p, 



(23) 



where p r := and where we recall that f(y) is the probability density of the jump 
Yi ~ N(^j,(Jj) : f(y) = (ctjv 7 ^-) -1 exp ( - (y - fij) 2 /2aj) in the Merton model, or 
/(y) = P«ie _aiX lx>o + (1 — p)oc2e° l2X 'ix<o in the Kou model. Before applying a suitable 



numerical integration algorithm to the integral terms in ( 23 ) , we truncate the integrals from 



an infinite computational range to a finite one. Briani et al (Briani et al., 2007), Cont and 



Voltchkova (Cont and Voltchkova, 2005), Tankov and Voltchkova (Tankov and Voltchkova 



2009) and d'Halluin et al. (d'Halluin et al., 2004 2005) have provided different numerical 



techniques to find out a finite computational range so as to reduce the numerical approxima- 
tion errors when doing this truncation. In this thesis we shall adopt the Briani et al. numerical 



technique to truncate the integral domain of our PIDE (cf. (Briani et al., 2007)) in both the 



Merton and Kou model. See [S] for a proof. Supposed e > 0, a formula of selecting a bounded 
interval [y~ e , y e ] for the set of points y in the Merton case is: 



V -""J 

V-e = ~Ve, Vy < 



2aj log(ecrjV2^/2) + fij, Vy > 



In the Kou model we have 



I/ e = log(e/p)/(l-«i), Vy >0 
y_ e = -log(e/(l-p))/(l-a 2 ) J Vy < 0, 



(24) 
(25) 



(26) 
(27) 



We therefore transform equation (23) into 

■2 



Ap T = y A xxP + ( r - q - y - Ary ) A, p + ( V + X)Ap + 
X^' A(y)f(y)dy ) j p. 



(28) 



We use matlab's adaptive Gauss-Kronrod quadrature to evaluate the matrix of the integrals 



in ([28]) : this amounts to approximating 

m 

4>{\xi + y - Xj\)f(y)dy « ^w fc </>(|xi + y k -Xj\)f(y k ), 



(29) 



k=l 



where and y^ are suitable quadrature weights and quadrature points; cf. (Shampine, 2008) 
for details. To simplify notations, we set 



m 

F(xi - Xj) = y^Wk^Xj + yk~ xj\)f(y k ). 

k=l 
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Then the integrals in equation (28) will be approximated by 



F(x 1 - xi) F(xi - x 2 ) 
F(x 2 - x\) F(x 2 - x 2 ) 



. F(xi - x N ) 
■ F(x 2 - x N ) 



F(x N - xi) F(x N -x 2 ) ... F(x N - x N )_ 
C(y). 



Substituting (30) into equation ( [28] ), we arrive at the new approximate equation: 

J2 



a' 



A xx p + r 



Xrj A x p + (r + X)Ap + XC(y)p. 



(30) 



(31) 



As we have known the Cubic Spline is strictly conditionally positive definite function of order 
2, the invertibility of A is not assumed without adding a real- valued polynomial of degree 



at most 1 in (16) (cf . | Wendland 2005). Nevertheless, Bos and Salkauskas proved that A is 
non-singular in a univariate case (cf. |Bos and Salkauskas 1987 Theorem 5.1). As a result, 
the invertibility of A is still guaranteed. 

Although the invertibility of A is able to be shown for all 4> of the interest, the inverse of A, 



may often be very ill-conditioned to solve when its size increases (cf. Fasshauer, 2007 



chapter 16). As a result, it may be impossible to solve accurately using standard floating 
point arithmetic. To address this problem, we factorise A into the following form (cf. Bos and 
Salkauskas[ |l987[ Theorem 3.7): 



A = FCF. 



(32) 



Here F is a N x N matrix, 



Fi 
\x 2 



Xl\ 



\Xl 

\x 2 



x 2 \ 
x 2 \ 



\Xl 

\x 2 



X3\ 
^3 



\x 2 



XN\ 
X N \ 



\XN — X\\ |Xjv — X 2 \ \XN — X3\ 

and C is a near tridiagonal N x N matrix, 



\XN ~ X N \ 



(33) 



7i 







q h 
2 

| 2h | 
| 2h | 







s 

2 







S 

2 








I 2h \ 
I h- 







s 



(34) 



where h is the distance between Xj+\ and Xj for 1 < i < N — 1 and S = Nh. We also have an 



explicit form of F (cf. Bos and Salkauskas 



1987 



Lemma 3.6) which is equal to 



h-S J_ 

_1_ 

2/f 



2h 




2h 







2h 
~h 







J_ 

25' 






1 

2S 



2h h 2h 
J_ h-S 
2h 2hS ■ 







(35) 



We perform Gaussian elimination with partial pivoting to calculate C 1 . Then, we multiply 
both sides of (31) by C 1 and F and we finally obtain the following homogeneous system 
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of ODEs with constant coefficients: 

/ 2 

p T =F- 1 C- 1 F- 1 (^-A xx +( 



r 



a 
~2 



\ V )A x + {r + X)A + XC(y))p 



Qp 



(36) 



where © is defined by the left hand side. After some numerical experimentation, we found 
that the matrix © is very stiff. To explain why © is stiff, we shall use the following example to 
illustrate it. Suppose we select our maximum and minimum logarithm price x mm ( log(S' m i n )) 
and x max (log(5 max )) in (22) equal to —10 and 10 respectively, then we use (22) to generate 



a list of 100 interpolation points. Based on the procedures and the ideas we have mentioned 



above we can get a 100 x 100 matrix © in (36). Then we measure the stiffness ratio of 0. The 



stiffness ratio is the quotient of the largest and the smallest eignvalues of the Jacobian matrix 
0. The ratio we have is 1.2864 x 10 5 . This implies that (36) is a stiff ODE and therefore 
we have to solve the ODEs by an implicit method, e.g. backward differentiation formulas 
(BDFs), a modified Rosenbrock formula of order 2, the trapezoidal rule or TR-BDF2, an 
implicit Runge-Kutta formula with a first stage that is a trapezoidal rule step and a second 
stage that is a backward differentiation formula of order two. In this paper we use former one. 



4. Numerical Results 



4-1 European Vanilla Options 



In this section we first present a simple scheme to construct our computational range. We 
then present the numerical results of our Cubic Spline approximation scheme and compare 
these with Black-Scholes, Merton and Kou's analytical option price formula for both puts and 
calls. Beside this, we also compare the results of our Cubic Spline approximation scheme with 
those of the Briani et al. finite difference method (FD) with implicit and explicit (IMEX) 



scheme in ( Briani et al. , 2007 ) and the Almendral et al. finite element method (FE) with 



backward differentiation formulas of order two (BDF2) and FD with BDF2 in (Almendral 
" 20051). 



and Oosterlee 



We use EMS (22) to choose our interpolation points. Based on this set of interpolation 
points, we can construct our computational range. We distribute the interpolation point uni- 
formly around the logarithm strike price, logK, in order to achieve a higher accuracy of 
pricing European Vanilla Option. Our scheme of distributing the interpolation point is shown 
in Figure [Tj The idea can be explained as follows: We set the range of [x mm , x max ] and then 
use EMS to create N interpolation points. We distribute the first N/2 points uniformly in 
[x min , log(AT)] and then the rest in [\og(K), x max ]. 



h 


I 


y 


X . 
min 


log (K) 


X 

max 



Figure 1. Uniform distributions of the interpolation points around the strike price by using EMS. The red dots are 
the interpolation points. The blue cross is the location of the logarithm strike price. 



In option trading the region of most interest is when the mean of the stock prices is close 
to the strike price. Typically, the probability for a stock to default or to be very far from the 
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strike price is small. Therefore we define the region of interest as follows: 

x t € [x min , x max ] := [log(K/20), log(2K) ]. (37) 

Based on this region, we can measure the accuracy of our RBF-approximation. We use a 
set of evaluation points xf" x , for which we will simply take the grid points 

Xi := xf- x = x miQ + jAx, j = 0, 1, 2, . . . , iV cval - 1. (38) 

Here Ax = (x max - x msx )/(N eval - 1) with x min < x min < x max < x max and N eval is the 
number of the evaluation points chosen. 

It is also of great interest to measure the rate of convergence of our Cubic Spline- 
approximation scheme. By defining At = 1/Mq, where Mq is the number of time steps and 
Ax = 1/N, where N is the number of interpolation points, we assume that 

E^ixuT) = C t (At) Rt + C x (Ax) R °° (39) 

for the max error and 

E 2 ( Xi , T) = C t {At) Rt + C x (Ax) R * (40) 

for the root-mean-square (rms) error. Here Eoc(xi,T) is the max error, E2(xi,T) is the rms 
error, Xi is i th evaluation point, T is the maturity time, both Ct and C% are constants, Rt is 
the rates of convergence in time and Roq and R2 are the rates of convergence in space. The 
formulae of calculating the max error and the rms errors are: 



E oo = 0< max \V{e x \r) -u(xi,T)\, (41) 



and 

Ei 



Y, \V(e^,T)-u(x h rW (42) 

^cval <i< Neml 



respectively. Here V(e x , r) and u(x, r) are the exact value and approximate value at the point 
(x, t) respectively. 

Since we compare the accuracy of our Cubic Spline-approximation scheme with that of 
FDM and FEM, we use the relative error, 



ttt /- \ \V(e x ,T) -u(x,t)\ 

Erel (X,T) = y (e a >T ) ' ( 43 ) 



as the measure of the accuracy. 



It is known (Merton, 1976) that the analytical price of a European call/put option in the 



Merton Jump-diffusion model is given by 

Vuj(S t ,r, K, r,q,a) 

-E eW " t f + " )T)k V B s(S,.r,K,r M) (44, 

where r = T — t is the time to maturity, rj = e^ J+ ^~ — 1 represents the expected percentage 
change in the stock price originating from a jump, a\ = a 2 + the observed volatility, 
rj, = r — \rj + k log(l + rj)/{T — t), q is the dividend and Vbs the Black-Scholes price of a call 
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and put, computed as 



VBs{St,r, K, r k ,a k ,< 



Ke- rkT $(d^ tk ) 



iTe- rfcT #(-d_ )fe ) - Ste- qT <S>(-d A 
where <&(•) is the cumulative normal distribution and 



-,k. 



call option, 
put option, 



,k 



log(S t /K)+(r k - 



g+^/2)T 



(L 



For the derivation of Vuj(St,T, K, r, q, a), we shall refer to the reader to (Merton, 1976, Cont 
and Tankovl [20041). 



In general, for models where the characteristic function of the Levy process is known, an 
analytical solution of PIDE ([£]) may be found using Fourier analysis (Carr and Madan, 1999 



Lewis, 2001). For the sake of simplicity and accuracy we propose Jackson et al.'s Fourier Space 



Time-Stepping method rather than Carr-Madan's Fast Fourier Transform (FFT) method 
(Carr and Madan, |1999 ) and Lewis's FFT method (Lewis, 2001). In brief, the idea of this 
method is based on the Fourier transform of the PIDE. By making use of FFT and inverse Fast 
Fourier transform (FFT -1 ), European Option price can be determined. The pricing formula 
of evaluating European option can be expressed as follows: 

V Kou (S, r,K,r, a, q) = FFT' 1 [FFT [V Kou (S,T)}e^ T ], (45) 

where ifj(z) is the characteristic function of the Kou model which can be defined as: 

/>u.i , (l—p)a 2 



a 2 z 2 



+ iz-f c + A( 



+ 



1 . 



2 x a\ — iz a>2 + iz 

and Vkou(S,T) is the payoff function ([l~0|). For more details of this method, we shall refer 



the reader to ( |Jackson et al. 2008). This method has been reported to have second order 
convergence in space in European cases. 



Our RBF-algorithm for numerically solving i9h with initial condition (10) runs as follows: 



'1) Find the RBF-approximation to the initial value u(x, 0) using ESM (see 22). This will 
provide us with a set of interpolation points x\, . . . , x n , together with an initial vector 
p(0) = (pi(0),...,piv(0)). 



(2) Then use p(0) as initial value for the system (36). By using any stiff ODE solver, we 
find out the p(T) at time T. 

(3) Finally, substitute p(T) back into J2f=i Pj (^)<MI X ~ x j\) to get an approximate value 
of u(x, T). 

In our numerical experiment we implement the algorithm in MATLAB R2007b. We select 
our maximum and minimum logarithm price x m \ n ( log(£ m i n )) and x max ( log(S , max )) , as be- 
fore, equal to —10 and 10 respectively. Because of achieving more accurate approximation of 

40 for finding a finite 



the integral in (28), we also set e in both 24 and 26 to be 3.72 x 10 

-e,Ue\- Moreover, we use function quadgk which implements adap- 



computational interval [y_ 
tive Gauss-Kronrod quadrature for computing equation (29) as well as function odel5s which 



implements backward differentiation formulas (BDFs) of order two for the calculation of equa- 



tion (36). The main reason of choosing it is the following: According to (Iserles 2009) BDFs 



of orders 1 and 2 are A-stable (the stability region includes the entire left half complex plane) 



Since (36) is stiff, according to Theorem 4.11 (The Dahlquist second barrier) of (Iserles, 2009), 
the highest order of an A stable multistep methocQ such as BDFs, is only two. We therefore 



1 Multistep methods are used for the numerical solution of ordinary differential equations. Conceptually, a numerical 
method starts from an initial point and then takes a short step forward in time to find the next solution point. The 
process continues with subsequent steps to map out the solution. 
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conclude that our solution is second order convergence in time. By setting Rt = 2, (39) and 



(40) become 

E^xuT) = C t (At) 2 + C x (Ax) R °° (46) 

E 2 (xi,T) = C t (At) 2 + C x (Ax) R * (47) 

for the European option. This conclusion is in line with the finding of ( Pettersson et al\ 
2008). In (Pettersson et al., 2008) Pettersson et al. show that second order in time can be 



achieved in a European case due to the second order time-stepping scheme, BDFs of order 2. 
Although they solve Black Schole PDE rather than PIDE in their paper, an similar approach of 
solving European option like our approximation scheme is applied. The rest of this section, we 
numerically show that and R 2 are equal to 2. Besides this, we will numerically approximate 
A and T and launch a comparison between our approximation scheme and FDM and FEM. 

All the parameters of all the tables except Table |3j [6] and [9] are chosen from different 
literature. The parameter a = 1 in Table [3j [6] and [9] is selected to stress our numerical 
algorithm. From Table [T] to [9j E^ and E 2 falls down when the number of the interpolation 
points N increases. Our Cubic Spline approximation scheme can get second order convergence 
in space. This is due to the limited smoothness of the Cubic Spline which has second order 
of convergence (cf. (Wendland, 2005)). In Figure [2] [3] and |4j oscillations do not occur around 
the strike K for small T when we approximate A and V. In Table 10, we compare the results 



of the FD used in Briani et al.'s paper (Briani et al., 2007) with those using our Cubic Spline 



approximation scheme. Our numerical approximation scheme can achieve lower E Te \_ (log S, T) 



than ARS-233 scheme and Explicit scheme. Table 10 and 12 are other comparisons of the 
accuracy between our Cubic Spline approximation scheme and Almendral and Oosterlee's 
FD and FE with BDF2. To illustrate a fair comparison, we set our maximum and minimum 
logarithm price x m i n and Almendral and Oosterlee proposed in their numerical 

experiments. Hence we set [x m i n rr max ] equal to [-4 4] and [-6 6] in the Merton model (Table 
11) and the Kou model (Table 12) respectively. Our Cubic Spline approximation scheme can 
attain lower E Te i (log S, T) than FD and FE with BDF2 in both the Merton and Kou cases. 



N 






E 2 (xi,T) 


R2 


100 


4.207101E-03 


N/A 


1.864736E-03 


N/A 


600 


1.195088E-04 


1.988 


5.143665E-05 


2.004 


1100 


3.554622E-05 


2.000 


1.528321E-05 


2.002 


1600 


1.679290E-05 


2.001 


7.219811E-06 


2.001 


2100 


9.745141E-06 


2.001 


4.189909E-06 


2.001 


2600 


6.354765E-06 


2.002 


2.732818E-06 


2.001 


3100 


4.468110E-06 


2.003 


1.921950E-06 


2.001 


3600 


3.311319E-06 


2.004 


1.424931E-06 


2.001 



Table 1. and E2 of the Cubic Spline approximation for pricing of a European put under the Black-Scholcs model are 

presented. iV is the number of the interpolation points. x% — log Si is any evaluation points of a range of S from 0.05 to 2 and 
the total numbers arc 1950. T is the Time-to-maturity. The parameters arc: r — 0.04, q — 0, tr = 0.29, K — 1 and T — 1. The 
parameters arc taken from Giles and Carter (2006). The order of convergence is 2 in space. 
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N 




Rod 




E 2 {x l ,T) 


R2 


100 


1.924131E-02 


N/A 


4. 


690135E-03 


N/A 


600 


7.143939E-04 


1.838 


1. 


.296858E-04 


2.003 


1100 


2.171519E-04 


1.965 


3. 


870772E-05 


1.995 


1600 


1.031950E-04 


1.986 


1. 


830673E-05 


1.998 


2100 


6.002721E-05 


1.992 


1. 


063352E-05 


1.998 


2600 


3.919766E-05 


1.995 


6. 


934013E-06 


2.002 


3100 


2.758717E-05 


1.997 


4. 


877540E-06 


2.000 


3600 


2.046213E-05 


1.998 


3. 


616699E-06 


2.000 



Table 2. Eac and E2 of the Cubic Spline approximation for pricing of a European call under the Black-Scholcs model are 
presented. iV is the number of the interpolation points. Xi — log Si is any evaluation points of a range of S from 0.05 to 2 and 
the total numbers are 1950. T is the Timc-to-maturity. The parameters arc: r — 0.05, q — 0, a — 0.2, K — 1 and T — 2. The 
parameters arc taken from Shaw (2009). The order of convergence is 2 in space. 



N 




Eocixii 1") 


Rod 




E 2 {x l ,T) 


R2 


100 


2 


.325676E-03 





.000 


1 


.40461 1E-03 


N/A 


600 


6 


.473617E-05 


1 


.999 


3 


.856043E-05 


2.007 


1100 


1 


.923322E-05 


2 


.002 


1 


.145625E-05 


2.002 


1600 


9 


.079037E-06 


2 


.003 


5 


.411921E-06 


2.001 


2100 


5 


.265272E-06 


2 


.004 


3 


.140776E-06 


2.001 


2600 


3 


.430306E-06 


2 


.006 


2 


.048580E-06 


2.001 


3100 


2 


.406208E-06 


2 


.016 


1 


.441039E-06 


2.000 


3600 


1 


.782442E-06 


2 


.007 


1 


.068202E-06 


2.002 



Table 3. E^, and E2 of the Cubic Spline approximation for pricing of a European call under the Black-Scholcs model are 
presented. iV is the number of the interpolation points, x% — log Si is any evaluation points of a range of S from 0.05 to 2 and 
the total numbers arc 1950. T is the Timc-to-maturity. The parameters are: r — 0.3, q — 0.1, a — 1. K — 1 and T — 0.25, 
whereas the parameter <r = 1 is selected to stress our numerical algorithm. The order of convergence is 2 in space. 



N 


Eoo(_Xi, T) 


Rco 


E 2 {xi,T) 


R2 


100 


1.428497E-02 


N/A 


3.749983E-03 


N/A 


600 


4.642130E-04 


1.912 


1.011341E-04 


2.016 


1100 


1.402519E-04 


1.975 


3.011378E-05 


1.999 


1600 


6.640377E-05 


1.995 


1.423346E-05 


2.000 


2100 


3.860331E-05 


1.995 


8.262241E-06 


2.000 


2600 


2.518672E-05 


1.999 


5.389115E-06 


2.001 


3100 


1.772559E-05 


1.997 


3.790660E-06 


2.000 


3600 


1.314288E-05 


2.000 


2.810697E-06 


2.000 



Table 4. E^, and E2 of the Cubic Spline approximation for pricing of a European call under the Mcrton Jump-diffusion model 
are presented. N is the number of the interpolation points. x% — log Si is any evaluation points of a range of S from 0.05 to 2 and 
the total numbers arc 1950. T is the Time-to-maturity. The parameters arc: r — 0.05, q — 0. a — 0.15, <j j — 0.45, fj,j — —0.9, 
A — 0.1, K — 1 and T — 0.25. The parameters are taken from ( And ersen and Andreasen| [2000[ |. The order of convergence is 2 
in space. 
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N 




Rod 


E 2 {x l ,T) 


R2 


100 


1.956920E-02 


N/A 


4.723349E-03 


N/A 


600 


7.326011E-04 


1.833 


1.305576E-04 


2.003 


1100 


2.240092E-04 


1.955 


3.898655E-05 


1.994 


1600 


1.069094E-04 


1.974 


1.844062E-05 


1.998 


2100 


6.223777E-05 


1.990 


1.071235E-05 


1.997 


2600 


4.062560E-05 


1.997 


6.985440E-06 


2.002 


3100 


2.859186E-05 


1.997 


4.913762E-06 


2.000 


3600 


2.121748E-05 


1.995 


3.643595E-06 


2.000 



Table 5. Kqo and E2 of the Cubic Spline approximation for pricing of a European put under the Mcrton Jump-diffusion model 
are presented. N is the number of the interpolation points. Xi — log Si is any evaluation points of a range of S from 0.05 to 
2 and the total numbers are 1950. T is the Timc-to-maturity. The parameters are: r — 0.05, q — 0.02, a — 0.15, aj — 0.4, 
fij — —1.08, A — 0.1, K — 1 and T — 0.1. The parameters are taken from | |Andersen and And rcascn 2000). The order of 
convergence is 2 in space. 



N 


Eoa(_Xi, 


Rod 


E 2 {x l ,T) 


R2 


100 


1.026524E-03 


N/A 


7.090253E-04 


N/A 


600 


2.819557E-05 


2.006 


1.945356E-05 


2.007 


1100 


8.415823E-06 


1.995 


5.762520E-06 


2.007 


1600 


3.999351E-06 


1.986 


2.712396E-06 


2.011 


2100 


2.373272E-06 


1.919 


1.559774E-06 


2.035 


2600 


1.601472E-06 


1.842 


1.004746E-06 


2.059 


3100 


1.136188E-06 


1.951 


7.021072E-07 


2.038 


3600 


8.358248E-07 


2.053 


5.221973E-07 


1.980 



Table 6. E^ and E2 of the Cubic Spline approximation for pricing of a European call under the Mcrton Jump-diffusion model 
arc presented. TV is the number of the interpolation points, afj — log Si is any evaluation points of a range of S from 0.05 to 2 and 
the total numbers arc 1950. T is the Time-to-maturity. The parameters arc: r — 0.05, q — 0.01, a — 1, uj — 0.6, fxj — —1.08, 
A — 0.1, K — 1 and T — 1, whereas the parameter a — 1 is selected to stress our numerical algorithm. The order of convergence 
is 2 in space. 



N 


Eooixii T) 


Rod 


E 2 {x l ,T) 


R2 


100 


1.239165E-02 


N/A 


3.422908E-03 


N/A 


600 


3.932126E-04 


1.926 


9.440247E-05 


2.004 


1100 


1.179555E-04 


1.986 


2.808850E-05 


2.000 


1600 


5.589111E-05 


1.993 


1.327392E-05 


2.000 


2100 


3.246588E-05 


1.998 


7.705266E-06 


2.000 


2600 


2.118103E-05 


2.000 


5.025765E-06 


2.001 


3100 


1.490021E-05 


2.000 


3.535171E-06 


2.000 


3600 


1.105067E-05 


1.999 


2.621377E-06 


2.000 



Table 7. Eqq and E2 of the Cubic Spline approximation for pricing of a European put under the Kou Jump-diffusion model 
are presented. TV is the number of the interpolation points. Xi — log is any evaluation points of a range of S from 0.05 to 2 
and the total numbers are 1950. T is the Time-to-maturity. The parameters are: r — 0, q — 0, <r — 0.2, a.± — 3, a>2 — 2, A — 0.2, 
p — 0.5, K — 1 and T — 0.2. The parameters arc taken from (Almcndral and Oostciicc 2005}. The order of convergence is 2 in 
space. 
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N 


Eoo(.Xii -0 


Rod 


E 2 {x l ,T) 


R2 


100 


1.433875E-02 


N/A 


3.766745E-03 


N/A 


600 


4.665677E-04 


1.912 


1.022079E-04 


2.013 


1100 


1.404381E-04 


1.981 


3.043034E-05 


1.999 


1600 


6.660275E-05 


1.991 


1.438190E-05 


2.000 


2100 


3.868283E-05 


1.998 


8.348098E-06 


2.000 


2600 


2.522395E-05 


2.002 


5.444331E-06 


2.001 


3100 


1.773247E-05 


2.003 


3.828943E-06 


2.001 


3600 


1.314079E-05 


2.004 


2.838628E-06 


2.001 



Table 8. Eqo an d &2 of the Cubic Spline approximation for pricing of a European call under the Kou Jump-diffusion model 
are presented. N is the number of the interpolation points. x% — log Si is any evaluation points of a range of S from 0.05 to 
2 and the total numbers are 1950. T is the Time-to-maturity. The parameters are: r — 0.05, q — 0, a — 0.15, ai — 3.0465, 
a 2 = 3.0465, A = 0.1, p = 0.3445, K = 1 and T = 0.25. The parameters are taken from jCarr and Mayo| [2007| |. The order of 
convergence is 2 in space. 



N 


Eoa{.Xi, 


Rod 


E 2 {x l ,T) 


R2 


100 


1.080306E-03 


N/A 


7.074108E-04 


N/A 


600 


2.973137E-05 


2.005 


1.940773E-05 


2.007 


1100 


8.870629E-06 


1.995 


5.757611E-06 


2.005 


1600 


4.229400E-06 


1.977 


2.712641E-06 


2.009 


2100 


2.490583E-06 


1.947 


1.567674E-06 


2.016 


2600 


1.674611E-06 


1.859 


1.014582E-06 


2.037 


3100 


1.191565E-06 


1.935 


7.096338E-07 


2.032 


3600 


9.018770E-07 


1.863 


5.232205E-07 


2.038 



Table 9. Eqq and E2 of the Cubic Spline approximation for pricing of a European put under the Kou Jump-diffusion model 
are presented. TV is the number of the interpolation points. £i — log Si is any evaluation points of a range of S from 0.05 to 2 
and the total numbers arc 1950. T is the Time-to-maturity. The parameters are: r — 0.04, q — 0.03, a — 1 , Qi — 4, 0.2 — 4, 
A — 0.3, p — 0.6 K — 1 and T — 1, whereas the parameter a — 1 is selected to stress our numerical algorithm. The order of 
convergence is 2 in space. 




0.5 1 1.5 2 0.5 1 1.5 2 



Figure 2. Put option Delta A (Left) and Gamma T (Right) in the Black-Scholes Model. The number of the interpolation 
points is 3600. The number of evaluation points of a range of S from 0.05 to 2 is 1950. The input parameters are provided 
in the caption to Table [3] 
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Figure 3. Call Option Delta A (Left) and Gamma V (Right)in the Merton Jump-diffusion Model. The number of the 
interpolation points is 3600. The number of evaluation points of a range of S from 0.05 to 2 is 1950. The input parameters 
are provided in the caption to Table [5] 




0.5 1 1.5 2 0.5 1 1.5 2 

S S 



Figure 4. Put Option Delta A (Left) and Gamma V (Right)in the Kou Jump-diffusion Model. The number of the 
interpolation points is 3600. The number of evaluation points of a range of S from 0.05 to 2 is 1950. The input parameters 
are provided in the caption to Table [8] 
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Explicit scheme ARS-233 Scheme 





N 


Value 


£ rel .(logS,T) 


Value 


E^log S,T) 


Call 


1024 


13.286915 


5.175624E-03 


13.287427 


5.214358E-03 


Put 


1024 


8.319940 


2.57797E-03 


8.326102 


1.839249E-03 






Cubic Spline 




N/A 




N 


Value 


£ rel .(logS,T) 


Value 


E vcl (logS,T) 


Call 


1024 


13.219358 


6.489263E-05 


N/A 


N/A 


Put 


1024 


8.342301 


1.027679E-04 


N/A 


N/A 



Table 10. Comparison between Explicit scheme i jBriani et aL | ( |2007) l ) , ARS-233 Scheme jBriani et aL|p007t )and Cubic Spline 
interpolation scheme in evaluating of European call/put under the Merton Jump-diffusion Model. The input parameters are: 
r = 0.05, q = 0, <t = 0.2, ctj = 0.8, /xj = 0, A = 0.1, K = 100, T = 1, and z = log 100. Reference prices of 13.218501 (call) and 
8.341444 (put) and parameters from |Briani et ar|j2007[ . 





FD with BDF2 


FE with BDF2 


N 
1025 


Value E Te i_(logS,T) 
9.411968E-02 1.682457e-04 


Value -E re i. (log S, T) 
9.412972E-02 6.165536E-05 




Cubic Spline 


N/A 


N 
1025 


Value E TeL (logS,T) 
9.413023E-02 5.621522E-005 


Value E Ie i (log S, T) 
N/A N/A 



Table 11. Comparison of FD with BDF2 jAlmendral and Oosterlee| | [2005l l), FE with BDF2 jAlmendral and Oosterleej J2005] >) 
and Cubic Spline interpolation scheme in evaluating of a European call (put) under the Merton Jump-diffusion Model. The 
input parameters are: r — 0, q — 0, er — 0.2, er,/ — 0.5, fij — 0, A — 0.1, K = 1, T = 1, and 3 = 1. Reference prices of 
0.094135525 for both call and put and parameters from [AlmendraT and Oostcrl ee] ( |2005*| > - 

FD with BDF2 FE with BDF2 



N 


Value E reL (logS,T) 


Value .Erei. (log S, T) 


513 


4.240E-02 6.346096E-03 


4.24579E-02 5.1285862E-03 




Cubic Spline 


N/A 


N 


Value E vc i(log S,T) 


Value E'rci. (log S, T) 


513 


4.254583E-02 3.061686E-03 


N/A N/A 



Table 12. Comparison of FD with BDF2 jAlmendral and Oosterlee| | [2005l l), FE with BDF2 jAlmendral and Oosterlee| p005l l) 
and Cubic Spline interpolation scheme in evaluating of a European call (put) under the Kou Jump-diffusion Model. The input 
parameters arc: r — 0, q — 0, cr — 0.2, cx\ — 3 : a 2 — 2, A — 0.2, p — 0.5, K — 1, T — 0.2, and S — 1. Reference prices of 
0.0426761 for both call and put and parameters from [Almendral and Oostcrlee (2005). 
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4-2 American Vanilla Put Options 

In this section we adapt an RBF-algorithm to compute American put-option prices. We then 
compare the option prices obtained from our RBF-algorithm with the Jackson et al. FST 
methods of (Jackson et al. 2008). As mentioned in Section [2j an American put option problem 



is a free boundary problem because of the possibility of early exercise at any point during its 
life, leading to the free boundary condition: 

u(x, t) = max (K — e x , u(x, r)j . 

Together with the smooth pasting condition mentioned in section [2j this uniquely determines 
the exercise boundary. 

The Jackson et al. FST methods suggest that their solutions can achieve second order in 
space when they implement their methods to price American put options. They implement 
their methods in the context of the LCP. As we have seen in Section[2j the value of an American 
option u(t,x) is always greater than or equal to the payoff function G(e x ). To numerically 
keep the condition u(t, x) — G(e x ) > to be continuously held (see Section [2]), this can be 
achieved when boundary conditions are applied. The numerical algorithm for this idea can be 
defined as follows: 

V(S, (m+l)At,K, r, a, q) 

= max{FFT _1 [ FFT [V(S, mAt, K, r, a, q) ]e M * ] , G(e x )} (48) 

where time interval At is obtained by dividing time-to-maturity T by the total number M, 
mA is the time-step, where m € {0, 1, 2, . . . , M — 1}, ip{z) is the characteristic function of the 
Merton/Kou models, V(S, (m + 1) At, K, r, a, q) is the American put price at time (m + 1) At 
and the payoff condition G(e x ) is equal to m&x(K — e x , 0). These methods also are required to 
swap between real and Fourier spaces at each time-step when the American option prices are 
calculated at each time interval. This is due to no convenient representation of the max(., .) 
operator in Fourier space. For the full schematic and numerical description of this method, 



we refer readers to (Jackson et al, 2008) 



As before, we use ESM to approximate u(x, 0) = ma.x(K — e x , 0) and then continue to work 
with the interpolation points found at r = 0. The algorithm now reads as follows: 

(1) Divide time-to-maturity T by total numbers of time-steps M to obtain time interval 
At and create a list of equally spaced time-points mAt, m £ {0, 1, 2, . . . , M — 1}. 

(2) Find the RBF-approximation to the initial value u(x, 0) using ESM. This will provide 
us with a set of interpolation points x\, .. . ,x n , together with an initial vector p(0) = 
( Pl (0),...,p N (0)) 



(3) Assume we have already determined p(mAt) (if m = 0, we have p(0)) in equation (36). 
Solve the system of (stiff) ODEs to find p((m + l)At) at the next successive time-step, 
(m + I) At. 

(4) Then at time {m + l)At, for each interpolation point Xi, define 

N 

u(xi, (m + l)At) = max ((K - e Xi ), ^Pj{{m + l)At)4>{\xi - xj\)). 

i=i 

(5) Find a new vector p((m+l) At) such that u(x{, (m+1) At) = J2f=i Pj {( m +l)A-t)(l)(\xi — 
Xj\) for all i. 

(6) Repeat Step 3.) to 5.) until m = M — 1. 

(7) Finally, substitute p(T) back into Ylf=i PjCr)(ft{\ x ~ x j\) to get an approximate value 



of u{x, T). 



The settings of our numerical experiment are the same as those in section 



4.1 



We also 
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(49) 



calculate the rate of convergence in time. If we hold Ax constant, (39) and (40) become 

£oc(x*,T) = C t (At) Rt 
and 

E 2 ( Xi ,T) = C t (At) R < (50) 

respectively. 



The results from Table 13 to 18 suggest that our Cubic Spline approximation method for 
pricing of American put options is second order in spatial variables and first order in time 
variables when the number of interpolation numbers N and the number of time-steps Mq 



are twofold and fourfold respectively. In table 19 and 20, we implement our own BDF2 with 



fixed time steps rather than using odelbs with variable time steps. From these two tables, 
we can achieve first order in time variables when the number of interpolation numbers N is 
held constant and the number of time-steps Mq is quadrupled. Moreover, from Figure [5] to [8j 
oscillations do not occur around the strike K for small or big T when we approximate A and 

r. 



N 


M 




Roo 


E 2 (xi,T) 


R 2 


225 


10 


2.368536E-03 


N/A 


1.007946E-03 


N/A 


450 


40 


7.746936E-04 


1.612 


2.740154E-04 


1.879 


900 


160 


2.260415E-04 


1.777 


6.969946E-05 


1.975 


1800 


640 


6.362341E-05 


1.829 


1.888980E-05 


1.884 


3600 


2560 


1.613907E-05 


1.979 


4.715908E-06 


2.002 



Table 13. E^ and E2 of the Cubic Spline approximation for pricing of an American put under the Mcrton model arc presented. 
N is the number of the interpolation points. Mo is the number of the time steps. Xi — log is any evaluation points of a range 
of S from 0.05 to 2 and the total numbers arc 1950. T is the Timc-to-maturity. The parameters arc: r — 0.05, q — 0, a — 0.15, 
crj — 0.45, f-i j — —0.9, A — 0.1, K — 1 and T — 0.25. The parameters arc taken from (Andersen and Andrcascn. 2000}. The 
order of convergence is 2 in space and 1 in time. 



N 


M 




Roo 


-^00 {%i ) T") 


R2 


225 


10 


3.401417E-03 


N/A 


7.995993E-04 


N/A 


450 


40 


1.318325E-03 


1.367 


2.451148E-04 


1.706 


900 


160 


3.744579E-04 


1.816 


6.873071E-05 


1.834 


1800 


640 


1.055849E-04 


1.826 


1.927219E-05 


1.834 


3600 


2560 


2.823205E-05 


1.903 


5.121082E-06 


1.912 



Table 14. E^ and E2 of the Cubic Spline approximation for pricing of an American put under the Mcrton model arc presented. 
N is the number of the interpolation points. Mo is the number of the time steps. Xi — log Si is any evaluation points of a range 
of S from 0.05 to 2 and the total numbers arc 1950. T is the Timc-to-maturity. The parameters arc: r — 0.05, q — 0.02, a — 0.15, 
ctj — 0.4, fjbj — —1.08, A — 0.1, K — 1 and T — 0.1. The parameters arc taken from (Andersen and Andrcascn. 2000). The 
order of convergence is 2 in space and 1 in time. 



5. Conclusion 



We have implemented an RBF method to solve the PIDE boundary value problem for pricing 
American put and European call/put options on a dividend-paying stock in the Merton and 
Kou Jump- diffusion market. By using the numerical scheme of Briani et ai., we find out a 
finite computational range of our global integral. Our results suggest that the Cubic Spline 
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N 


M 


Eoo{Xii 


Roo 


E 2 (x u T) 


R2 


225 


10 


4.935878E-03 


N/A 


1.613323E-03 


N/A 


450 


40 


1.236617E-03 


1.997 


3.725615E-04 


2.114 


900 


160 


3.093198E-04 


1.999 


9.101657E-05 


2.033 


1800 


640 


7.734030E-05 


2.000 


2.133679E-05 


2.093 


3600 


2560 


1.932168E-005 


2.001 


5.074520E-06 


2.072 



Table 15. E^ and E2 of the Cubic Spline approximation for pricing of an American put under the Mcrton model arc presented. 
N is the number of the interpolation points. M$ is the number of the time steps, x^ — log Si is any evaluation points of a range 
of S from 0.05 to 2 and the total numbers arc 1950. T is the Timc-to-maturity. The parameters arc: r — 0.05, q — 0.01, <T = 1, 
ctj — 0.6, pj — —1.08, A — 0.1, K — 1 and T — 1, whereas the parameter a = 1 is selected to stress our numerical algorithm. 
The order of convergence is 2 in space and 1 in time. 



N 


M 


Eoo (,X{, T") 


Roo 


E 2 (xi,T) 


R2 


225 


10 


1.508321E-03 


N/A 


5.589125E-04 


N/A 


450 


40 


7.233939E-04 


1.060 


1.759571E-04 


1.667 


900 


160 


1.958968E-04 


1.885 


4.733738E-05 


1.894 


1800 


640 


5.243753E-05 


1.901 


1.271703E-05 


1.896 


3600 


2560 


1.374207E-05 


1.932 


3.405083E-06 


1.901 



Table 16. E^ and E2 of the Cubic Spline approximation for pricing of an American put under the Kou Jump-diffusion model 
are presented. TV is the number of the interpolation points. Mo is the number of the time steps. Xi — log is any evaluation 
points of a range of S from 0.05 to 2 and the total numbers arc 1950. T is the Timc-to-maturity. The parameters are: r — 0, 
q — 0, <7 — 0.2, Qi — 3, a>2 — 2, A — 0.2, p — 0.5, K — 1 and T — 0.2. The parameters arc taken from ( Almcndral and Oostcrlcc 
|2005[|. The order of convergence is 2 in space and 1 in time. 



N 


M 


Eoo {Xi, T") 


Roo 


E 2 (xi,T) 


R2 


225 


10 


1.933354E-03 


N/A 


8.983577E-04 


N/A 


450 


40 


8.487095E-04 


1.188 


2.783005E-04 


1.691 


900 


160 


2.497213E-04 


1.765 


7.257535E-05 


1.939 


1800 


640 


6.843085E-05 


1.868 


1.933309E-05 


1.908 


3600 


2560 


1.827216E-05 


1.905 


5.119491E-06 


1.917 



Table 17. E^ and E2 of the Cubic Spline approximation for pricing of an American put under the Kou Jump-diffusion model 
are presented. TV is the number of the interpolation points. Mo is the number of the time steps. a% — log is any evaluation 
points of a range of S from 0.05 to 2 and the total numbers are 1950. T is the Timc-to-maturity. The parameters are: r — 0.05, 
q = 0, a = 0.15, at = 3.0465, ot 2 = 3.0465, A = 0.1, p = 0.3445, K — 1 and T = 0.25. The parameters are taken from | |Carr| 
|and Mayo"l|2007) . The order of convergence is 2 in space and 1 in time. 



N 


M 


Eoo (.Xi, T") 


Roo 


E 2 (xi,T) 


R2 


225 


10 


3.839148E-03 


N/A 


1.095217E-03 


N/A 


450 


40 


9.616353E-04 


1.997 


2.458977E-04 


2.155 


900 


160 


2.405238E-04 


1.999 


6.111403E-05 


2.008 


1800 


640 


6.013812E-05 


2.000 


1.508359E-05 


2.019 


3600 


2560 


1.490999E-05 


2.012 


3.768285E-06 


2.001 



Table 18. E^ and E2 of the Cubic Spline approximation for pricing of an American put under the Kou Jump-diffusion model 
are presented. TV is the number of the interpolation points. Mo is the number of the time steps. Xi — log Si is any evaluation 
points of a range of S from 0.05 to 2 and the total numbers are 1950. T is the Timc-to-maturity. The parameters are: r — 0.04, 
q — 0.03, (7 — 1, Q'i — 4, a 2 — 4, A — 0.3, p — 0.6, K — 1 and T — 1, whereas the parameter a — 1 is selected to stress our 
numerical algorithm. The order of convergence is 2 in space and 1 in time. 
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N 


M 






E 2 (xi,T) 


R2 


3600 


A C\ 

40 


1.029993E-03 


TvT / A 

JN /A 


3.556438E-04 


AT / A 

JN /A 


3600 


160 


2.974325E-04 


0.896 


8.866477E-05 


1.002 


3600 


640 


8.273457E-05 


0.923 


2.147058E-05 


1.023 


3600 


2560 


2.123568E-05 


0.982 


5.345367E-06 


1.003 


3600 


40 


1.523857E-03 


N/A 


3.755414E-04 


N/A 


3600 


160 


4.561966E-04 


0.870 


1.108779E-04 


0.880 


3600 


640 


1.253234E-04 


0.932 


2.857790E-05 


0.978 


3600 


2560 


3.230111E-05 


0.978 


7.134578E-06 


1.001 



Table 19. E^o and E2 of the Cubic Spline approximation for pricing of an American put under the Mcrton Jump-diffusion 
model are presented. TV is the number of the interpolation points. Mo is the number of the time steps. Xi — log Si is any 
evaluation points of a range of 5 from 0.05 to 2 and the total numbers arc 1950. Top: The input parameters are provided in the 
caption to Tabic [l3| Bottom: The input parameters are provided in the caption to Table [l4| The order of convergence is 1 in 
time. 



N 


M 


Eoo (-^i) -0 




E 2 (xi,T) 


R2 


3600 


40 


1.687267E-03 


N/A 


1.878158E-04 


N/A 


3600 


160 


4.913041E-04 


0.890 


5.514562E-05 


0.884 


3600 


640 


1.366624E-04 


0.923 


1.590242E-05 


0.897 


3600 


2560 


3.478690E-05 


0.987 


3.768285E-06 


0.923 


3600 


40 


8.987023E-04 


N/A 


2.671173E-04 


N/A 


3600 


160 


3.151073E-04 


0.756 


7.544376E-05 


0.912 


3600 


640 


9.355218E-05 


0.876 


2.184649E-05 


0.894 


3600 


2560 


2.675431E-05 


0.903 


6.204563E-06 


0.908 



Table 20. and E2 of the Cubic Spline approximation for pricing of an American put under the Kou Jump-diffusion model 

are presented. TV is the number of the interpolation points. Mo is the number of the time steps. a% — log Si is any evaluation 
points of a range of S from 0.05 to 2 and the total numbers arc 1950. Top: The input parameters arc provided in the caption to 
Tabic [Trjj Bottom: The input parameters are provided in the caption to Tabic [l7| The order of convergence is 1 in time. 



< 




0.5 1 1.5 2 0.5 1 1.5 2 

S S 



Figure 5. Put option Delta A (Left) and Gamma T (Right) in the Merton Jump diffusion Model. The number of the 
interpolation points N is 1800 and the number of time steps Mo is 640. The numb er o f evaluation points of a range of 
S from 0.05 to 2 is 1950. The input parameters are provided in the caption to Table fl3] 
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Figure 6. Put option Delta A (Left) and Gamma T (Right) in the Merton Jump diffusion Model. The number of the 
interpolation points N is 1800 and the number of time steps Mq is 640. The number of evaluation points of a range of 
S from 0.05 to 2 is 1950. The input parameters are provided in the caption to Table [l4| 





Figure 7. Put Option Delta A (Left) and Gamma T (Right)in the Kou Jump-diffusion Model. The number of the 
interpolation points N is 1800 and the number of time steps Mo is 640. The number of evaluation points of a range of 
S from 0.05 to 2 is 1950. The input parameters are provided in the caption to Table [l7| 



approximation scheme can achieve second-order convergence in both spatial and time vari- 
ables (due to the second order time-stepping scheme, BDFs of order 2) when it is used to 
compute European call/put options. Moreover, the results also show that our approximation 
solution can get second-order convergence in spatial variables and first-order convergence in 
time variables when the approximation scheme is used to compute American put options. 
Beside this, we compare our RBF-approximation method against FDM and FEM. Our results 
suggest that one can achieve a high accuracy by implementing our meshless scheme. Moreover, 
in terms of meshless interpolation methods, we use cubic spline as a basis function rather than 
MQ. This basis function can avoid the open question of choosing an optimal shape parameter 
of MQ. Beside this, by using factorisation of the Cubic Spline, we can avoid inverting an ill- 
conditioned cubic spline interpolant directly. Finally, throughout the analysis of both A and 
r, our RBF-approximation method can also avoid the oscillation problem around the strike 
K in both American and European cases. 
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0.5 1 1.5 2 0.5 1 1.5 2 

S S 



Figure 8. Put Option Delta A (Left) and Gamma V (Right)in the Kou Jump-diffusion Model. The number of the 
interpolation points N is 1800 and the number of time steps Mo is 640. The number of evaluation points of a range of 
S from 0.05 to 2 is 1950. The input parameters are provided in the caption to Table [l8| 

At this stage of development, the Cubic Spline approximation scheme is first order in time 
for American put options although a second order time-stepping scheme, BDFs of order 2 is 
implemented. We are investigating various approaches to improve the Cubic Spline approx- 
imation for time variables and will treat them in a future paper. Our Method extends in 
principle to pure jump Levy type models for the underlying stocks, like the Variance Gamma 
(VG) model or the CGMY model. 
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Appendix A. A Finite Computational Range in the Jump-diffusion Model 

In the Merton Model suppose in a domain f2 E K European option price u(x, r) satisfies 
Lipchitz inequality such that 

\u(xi, t) — u(x2, t)| < L\x\ — X%\, \/x\,X2 G O. 

Then we choose a parameter e > and select the bounded intervals , y e ] as the set of all 
points y that verify 



k(y) 



2naj 



'j > e. 



Because of the symmetry of k(y) we set y_ e = — y e . Then the truncation of the integral domain 
giving an error to approximation of the problem can be estimated by 



{u{x + y) - u{x))k{y) dy - / {u(x + y) - u{x))k{y) dy 



< L 



< L 



(x + y - x)k(y) dy - j (x + y - x)k{y) dy 



00 



fOO 



y\k{y)dy+ I \y\k{y)dy 



1 



exp( 



2aj 



'dy 



f°° 1 y 

2 / (y + fJ-j) nr- exp(-- T )dy 



1 y z 
2 / (y + fJ-j) exp(-- T )dy 
'v.-nj V2iraj 2a J 



< 2 



1 



y 



2/ e -AO 



(2/ + y)—f= — ex p(-^2 ) d y 



2i\a j 



4trj (y e - /ij) 2 , 

: exp(- 



2vr 

2^6 



2aj 



(Ala) 

(Alb) 

(Ale) 

(Aid) 

(Ale) 

(Alf) 

(Alg) 
(Alh) 



Hence by using (Alg) and (Alh), 



//« - \l -2aj log(ecTjV27r/2) + 



(A2) 



We use the aforementioned arguments to find the finite computational range [y— e , y e ] m the 
Kou model. We carry out the reasoning for the positive semi-axis (the reasoning goes similarly 
for the negative semi-axis) and set k(y) = pa\e~ aiV for y > ((1 — p)a2e a2X for y < 0). Then, 
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y £ can be found out by the following equations: 



< L 



oo ry t 

(u(x + y)-u(x))Xf(y)dy- / (u(x + y) - u(x))Xf(y) dy 
o Jo 

oo ry e 

(x + y - x)Xf(y) dy- (x + y - x)Xf(y) dy 
o J o 



poo 

<L \y\f(y)dy 



\y\pot\e aiV dy 



pa\e 



(Gradshteyn and Ryzhik, 1994, equation 3.351 



V_ 
a i 



3/e(l-0!l) 



= pe 
= e, 

as a result, 

y e = log(e/p)/(l - ai). 
Similar arguments can be applied to y < 0, so 

y- e = -log (e/(l-p))/(l-a 2 ). 



(A3a) 
(A3b) 
(A3c) 
(A3d) 

(A3e) 

(A3f) 

(A3g) 
(A3h) 

(A4) 

(A5) 
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